knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE,
fig.height = 4, fig.width = 7)
library(terra)
library(sf)
library(data.table)
library(tidyverse)
library(here)
library(oharac)
source(here('common_fxns.R'))
Using functional entities defined in a prior script based on a suite of categorical traits, in combination with spatial distribution of species from AquaMaps and IUCN, calculate functional diversity metrics based on Mouillot et al (2014) and code from Sebastien Villeger: functional richness, functional vulnerability, functional redundancy, and functional over-redundancy.
Included species are:
spp_info <- assemble_spp_info_df() %>%
select(species, taxon, fe_id, mob, trp, len, wcol) %>%
distinct() %>%
mutate(mob = factor(mob, levels = c('ses', 'sed', 'mob', 'mig')),
trp = factor(trp, levels = c('pri', 'sec', 'ter', 'apex')),
len = factor(len, levels = c('tiny', 'vsm', 'sm', 'med', 'ml', 'lg', 'vlg', 'huge')))
n_fes_possible <- n_distinct(spp_info$mob) *
n_distinct(spp_info$trp) *
n_distinct(spp_info$len) *
n_distinct(spp_info$wcol)
n_fes_realized <- spp_info$fe_id %>% n_distinct()
Overall, with this set of traits used to define functional entities, we could have up to 512 FEs, only 339 (66.2%) of which are realized in this set of species. Note that the Mouillot 2014 paper includes two additional traits bringing the total number of possible FEs up to 5670, only 646 of which are realized by the 6316 tropical fish spp in their study.
knitr::kable(t(table(spp_info$len)))
| tiny | vsm | sm | med | ml | lg | vlg | huge |
|---|---|---|---|---|---|---|---|
| 2174 | 5959 | 4179 | 3595 | 2555 | 1370 | 906 | 406 |
knitr::kable(t(table(spp_info$trp)))
| pri | sec | ter | apex |
|---|---|---|---|
| 1647 | 4247 | 12810 | 2440 |
knitr::kable(t(table(spp_info$mob)))
| ses | sed | mob | mig |
|---|---|---|---|
| 2539 | 4378 | 10709 | 3518 |
knitr::kable(t(table(spp_info$wcol)))
| ben | bp | pel | rf |
|---|---|---|---|
| 12702 | 1130 | 2233 | 5079 |
Here we use the IUCN species range maps rasterized to 10 km
Mollweide, and AquaMaps species distributions reprojected to 10 km
Mollweide. Priority: IUCN species maps, then AquaMaps with occur_cells ≥
10. These map filenames are already set up in the
assemble_spp_info_df() function.
Read in functional entity assignments and join to spatial ranges. Per cell per FE, calculate the number of species for later calculation into functional richness, redundancy, vulnerability. Break the process into chunks to avoid overloading memory.
process_spp_to_fe <- function(fe, spp_maps) { ### fe <- fe_vec[1]
dt_out <- spp_maps %>%
filter(fe_id == fe) %>%
data.table() %>%
.[ , .(n_spp_fe = length(unique(species))),
by = .(cell_id, fe_id)]
return(dt_out)
}
fe_summary_file <- here_anx('func_entities', 'fe_cell_summary.csv')
# unlink(fe_summary_file)
if(!file.exists(fe_summary_file)) {
spp_info_df <- assemble_spp_info_df(fe_only = TRUE, vuln_only = TRUE)
### table(spp_info_df %>% select(species, len) %>% distinct() %>% .$len)
# huge lg med ml sm tiny vlg vsm
# 406 1370 3595 2555 4179 2174 906 5959
### Chunk out and save smaller files to avoid loading in EVERY spp map
### all at once!
n_chunks <- 40
spp_to_chunks <- spp_info_df %>%
select(species) %>%
distinct() %>%
mutate(chunk = rep(1:n_chunks, length.out = n()))
spp_clean <- spp_info_df %>%
select(species, fe_id, map_f, taxon) %>%
distinct() %>%
left_join(spp_to_chunks, by = 'species')
spp_fes <- spp_clean %>%
select(species, fe_id) %>%
distinct()
tmp_chunk_filestem <- here('tmp/fe_summary_chunk_%s.csv')
for(i in 1:n_chunks) { ### i <- 1
tmp_chunk_f <- sprintf(tmp_chunk_filestem, i)
if(file.exists(tmp_chunk_f)) {
message('Temp file ', basename(tmp_chunk_f), ' exists... skipping!')
next()
}
message('Assembling species maps and binding FEs for chunk ', i, ' of ', n_chunks, '...')
spp_chunk <- spp_clean %>% filter(chunk == i)
chunk_spp_maps <- collect_spp_rangemaps(spp_vec = spp_chunk$species,
file_vec = spp_chunk$map_f,
idcol = 'species',
) %>%
dt_join(spp_fes, by = 'species', type = 'left') %>%
distinct()
fe_vec <- chunk_spp_maps$fe_id %>% unique() %>% sort()
message('... summarizing FR map for chunk ', i, ': ', length(fe_vec), ' FEs in ',
nrow(chunk_spp_maps), ' cell observations...')
chunk_fe_sumlist <- parallel::mclapply(fe_vec, mc.cores = 20,
FUN = process_spp_to_fe,
spp_maps = chunk_spp_maps)
if(check_tryerror(chunk_fe_sumlist)) {
stop('Try-errors detected in FE summary loop!')
}
message('... binding chunk and writing chunk summary to ', basename(tmp_chunk_f), '...')
chunk_fe_summary <- chunk_fe_sumlist %>% data.table::rbindlist()
fwrite(chunk_fe_summary, tmp_chunk_f)
}
message('... reading temp chunks, binding, and summarizing by fe and cell...')
fe_tmp_fs <- list.files(dirname(tmp_chunk_filestem),
pattern = 'fe_summary_chunk', full.names = TRUE)
fe_summary_list <- parallel::mclapply(fe_tmp_fs, FUN = fread, mc.cores = 10)
fe_summary <- rbindlist(fe_summary_list) %>%
.[ , .(n_spp_fe = sum(n_spp_fe)),
by = .(cell_id, fe_id)]
message('... writing summary to ', basename(fe_summary_file), '...')
fwrite(fe_summary, fe_summary_file)
### unlink(fe_tmp_fs)
}
Calculate functional diversity metrics from IUCN and AquaMaps summaries of functional entity membership per cell.
This set of calculations replicates the results of Villeger’s function but without resorting to a presence/absence matrix, instead just relying on the original AquaMaps species-cell table.
fe_metrics_file <- here_anx('func_entities/fe_metrics_per_cell.csv')
### unlink(fe_metrics_file)
if(!file.exists(fe_metrics_file)) {
# system.time({
fe_sum_total <- data.table::fread(fe_summary_file)
### there are 6.5 million cells. Chop into 500k instances and mclapply to summarize
# cell_id_vec <- 1:ncell(raster::raster(here('_spatial/ocean_area_mol.tif')))
chunk_size <- 250000
n_chunks <- 6.5e6 / chunk_size
nspp_tmp_f <- here('tmp/nspp_df_tmp.csv')
if(!file.exists(nspp_tmp_f)) {
message('Calculating species richness map...')
nspp_list <- parallel::mclapply(1:n_chunks, mc.cores = ceiling(n_chunks / 4),
FUN = function(n) { ### n <- 1
message('... spp richness chunk ', n, ' of ', n_chunks, '...')
cell_id_min <- (n - 1) * chunk_size + 1
cell_id_max <- n * chunk_size
nspp_sum <- fe_sum_total %>%
filter(between(cell_id, cell_id_min, cell_id_max)) %>%
data.table() %>%
.[ , .(n_spp = sum(n_spp_fe, na.rm = TRUE),
n_fe = n_distinct(fe_id)),
by = 'cell_id']
})
if(check_tryerror(nspp_list)) stop('Try errors detected!')
nspp_df <- nspp_list %>%
data.table::rbindlist()
write_csv(nspp_df, nspp_tmp_f)
} else {
message('File exists: ', basename(nspp_tmp_f))
nspp_df <- data.table::fread(nspp_tmp_f)
}
fvuln_tmp_f <- here('tmp/fvuln_df_tmp.csv')
if(!file.exists(fvuln_tmp_f)) {
message('Calculating functional vulnerability map...')
fv_list <- parallel::mclapply(1:n_chunks, mc.cores = ceiling(n_chunks / 4),
FUN = function(n) { ### n <- 1
message('... funct vuln chunk ', n, ' of ', n_chunks, '...')
cell_id_min <- (n - 1) * chunk_size + 1
cell_id_max <- n * chunk_size
f_vuln <- fe_sum_total %>%
filter(between(cell_id, cell_id_min, cell_id_max)) %>%
oharac::dt_join(nspp_df, by = 'cell_id', type = 'left') %>%
data.table() %>%
.[ , .(f_vuln = sum(n_spp_fe == 1) / n_fe,
f_wvuln = mean(0.5^(n_spp_fe - 1)),
f_red = mean(n_spp_fe)),
by = .(cell_id, n_fe, n_spp)]
})
if(check_tryerror(fv_list)) stop('Try errors detected!')
f_vuln_df <- fv_list %>%
data.table::rbindlist()
write_csv(f_vuln_df, fvuln_tmp_f)
} else {
message('File exists: ', basename(fvuln_tmp_f))
f_vuln_df <- data.table::fread(fvuln_tmp_f)
}
fred_tmp_f <- here('tmp/f_red_df_tmp.csv')
if(!file.exists(fred_tmp_f)) {
message('Calculating functional overredundancy map...')
f_overred_list <- parallel::mclapply(1:n_chunks, mc.cores = ceiling(n_chunks / 4),
FUN = function(n) { ### n <- 1
message('... funct overredundancy chunk ', n, ' of ', n_chunks, '...')
cell_id_min <- (n - 1) * chunk_size + 1
cell_id_max <- n * chunk_size
f_overred <- fe_sum_total %>%
filter(between(cell_id, cell_id_min, cell_id_max)) %>%
oharac::dt_join(f_vuln_df, by = 'cell_id', type = 'left') %>%
data.table() %>%
.[ , f_over := ifelse(n_spp_fe > f_red, n_spp_fe - f_red, 0)] %>%
.[ , .(f_overred = sum(f_over) / sum(n_spp_fe)),
by = 'cell_id']
})
if(check_tryerror(f_overred_list)) stop('Try errors detected!')
f_overred_df <- f_overred_list %>%
data.table::rbindlist()
write_csv(f_overred_df, fred_tmp_f)
} else {
message('File exists: ', basename(fred_tmp_f))
f_overred_df <- data.table::fread(fred_tmp_f)
}
message('Combining functional diversity metrics maps...')
fe_metrics_df <- f_vuln_df %>%
oharac::dt_join(f_overred_df, by = 'cell_id', type = 'full')
message('Combined set: ', nrow(fe_metrics_df), ' rows!')
write_csv(fe_metrics_df, fe_metrics_file)
}
Create and save rasters of functional diversity metrics
fe_metrics_df <- data.table::fread(fe_metrics_file)
### note that while there are 3,684,273 cells in this dataframe,
### 4273 of those are the Caspian Sea, so when those are dropped,
### it results in 3,680,000 exactly... ugh, weird
n_spp_rast <- map_to_mol(fe_metrics_df, which = 'n_spp')
n_fe_rast <- map_to_mol(fe_metrics_df, which = 'n_fe')
f_vuln_rast <- map_to_mol(fe_metrics_df, which = 'f_vuln')
f_wvuln_rast <- map_to_mol(fe_metrics_df, which = 'f_wvuln')
f_red_rast <- map_to_mol(fe_metrics_df, which = 'f_red')
f_overred_rast <- map_to_mol(fe_metrics_df, which = 'f_overred')
writeRaster(n_spp_rast, here('_output/func_entities/n_spp.tif'),
overwrite = TRUE)
writeRaster(n_fe_rast, here('_output/func_entities/n_fe.tif'),
overwrite = TRUE)
writeRaster(f_vuln_rast, here('_output/func_entities/f_vuln.tif'),
overwrite = TRUE)
writeRaster(f_wvuln_rast, here('_output/func_entities/f_wvuln.tif'),
overwrite = TRUE)
writeRaster(f_red_rast, here('_output/func_entities/f_red.tif'),
overwrite = TRUE)
writeRaster(f_overred_rast, here('_output/func_entities/f_overred.tif'),
overwrite = TRUE)
Examine distributions of each raster.
Plot each raster. NOTE: These no longer include spp that are missing
from vulnerability calculations - so the maps in
_output/nspp_maps (vuln \(\cap\) FE) should be identical to these in
output/func_entities.